{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 1. 以Lena为原始图像，通过OpenCV实现平均滤波，高斯滤波及中值滤波，比较滤波结果。\n",
    "#### (1) 原图与不同窗口大小的平均滤波的比较 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "def my_subplot(rows, cols, seq, img, title):\n",
    "    plt.subplot(rows, cols, seq)\n",
    "    plt.imshow(reverse(img))\n",
    "    plt.title(title)\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "\n",
    "\n",
    "src_img = cv.imread(r\"image/lena.jpg\")\n",
    "plt.figure(figsize = (2.5 * src_img.shape[0] / 100, 3.5 * src_img.shape[1] / 100))\n",
    "\n",
    "my_subplot(3, 2, 1, src_img, \"Original\")\n",
    "\n",
    "# 实现平均滤波 -- cv.blur(3*3)\n",
    "averaging_blur_result_1 = cv.blur(src_img, (3, 3))\n",
    "my_subplot(3, 2, 3, averaging_blur_result_1, \"Averaging Blur cv.blur(3*3)\")\n",
    "\n",
    "# 实现平均滤波 -- cv.boxFilter(3*3)\n",
    "averaging_blur_result_2 = cv.boxFilter(src_img, -1, (3, 3))\n",
    "my_subplot(3, 2, 4, averaging_blur_result_2, \"Averaging Blur cv.boxFilter(3*3)\")\n",
    "\n",
    "# 实现平均滤波 -- cv.blur(7*7)\n",
    "averaging_blur_result_3 = cv.blur(src_img, (7, 7))\n",
    "my_subplot(3, 2, 5, averaging_blur_result_3, \"Averaging Blur cv.blur(7*7)\")\n",
    "\n",
    "# 实现平均滤波 -- cv.filter2D(7*7)\n",
    "kernel = np.ones((7, 7), np.float32) / 49\n",
    "averaging_blur_result_4 = cv.filter2D(src_img, -1, kernel)\n",
    "my_subplot(3, 2, 6, averaging_blur_result_4, \"Averaging Blur cv.filter2D(7*7)\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果：\n",
    "![原图和平均滤波的对比](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/original_averaging_blur.png)\n",
    "#### 结果分析:\n",
    "+ Lena图像经过滤波后，灰度变化更加平滑；\n",
    "+ 滤波窗口的尺寸越大，Lena图像的滤波结果越模糊；"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (2) 原图与不同窗口大小的高斯滤波结果的比较"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "def my_plot(rows, cols, seq, img, title):\n",
    "    plt.subplot(rows, cols, seq)\n",
    "    plt.title(title)\n",
    "    plt.imshow(reverse(img))\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "\n",
    "    \n",
    "src_img = cv.imread(r\"image/lena.jpg\")\n",
    "\n",
    "plt.figure(figsize = (2.5 * src_img.shape[0] / 100, 3.5 * src_img.shape[1] / 100))\n",
    "\n",
    "# 画出原图\n",
    "my_plot(3, 2, 1, src_img, \"Original\")\n",
    "\n",
    "# 高斯窗口设置为3*3时\n",
    "gaussian_blur_result_3_3 = cv.GaussianBlur(src_img, (3, 3), sigmaX=0, sigmaY=0)\n",
    "my_plot(3, 2, 3, gaussian_blur_result_3_3, \"Gaussian Blur(3*3)\")\n",
    "\n",
    "# 高斯窗口设置为5*5时\n",
    "gaussian_blur_result_5_5 = cv.GaussianBlur(src_img, (5, 5), sigmaX=0, sigmaY=0)\n",
    "my_plot(3, 2, 4, gaussian_blur_result_5_5, \"Gaussian Blur(5*5)\")\n",
    "\n",
    "# 高斯窗口设置为7*7时\n",
    "gaussian_blur_result_7_7 = cv.GaussianBlur(src_img, (7, 7), sigmaX=0, sigmaY=0)\n",
    "my_plot(3, 2, 5, gaussian_blur_result_7_7, \"Gaussian Blur(7*7)\")\n",
    "\n",
    "# 高斯窗口设置为31*31时\n",
    "gaussian_blur_result_31_31 = cv.GaussianBlur(src_img, (31, 31), sigmaX=0, sigmaY=0)\n",
    "my_plot(3, 2, 6, gaussian_blur_result_31_31, \"Gaussian Blur(31*31)\")\n",
    "\n",
    "plt.show()\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果:\n",
    "![高斯滤波结果与原图的比较图](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/original_gaussian_blur.png)\n",
    "#### 结果分析:\n",
    "+ 将3\\*3高斯窗口滤波结果与原图比较，灰度过渡更为平滑，而且基本保留了原图的细节，效果较好；\n",
    "+ 将3\\*3，5\\*5，7\\*7，31\\*31高斯窗口滤波结果，与原图比较，可以看到，随着滤波窗口的变大，虽然灰度过渡越来越平滑，但窗口较大时，滤波结果显得越来越模糊，甚至丢失了原图的细节，可见窗口尺寸的选择不宜过大。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (3) 原图与不同窗口大小的中值滤波结果的比较"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "def my_plot(rows, cols, seq, img, title):\n",
    "    plt.subplot(rows, cols, seq)\n",
    "    plt.title(title)\n",
    "    plt.imshow(reverse(img))\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "\n",
    "    \n",
    "src_img = cv.imread(r\"image/lena.jpg\")\n",
    "plt.figure(figsize = (2.5 * src_img.shape[0] / 100, 2.5 * src_img.shape[1] / 100))\n",
    "\n",
    "# 画出原图\n",
    "my_plot(2, 2, 1, src_img, \"Original\")\n",
    "\n",
    "# 画出3*3中值滤波的结果\n",
    "median_blur_result_3_3 = cv.medianBlur(src_img, 3)\n",
    "my_plot(2, 2, 3, median_blur_result_3_3, \"Median Blur(3*3)\")\n",
    "\n",
    "# 画出5*5中值滤波的结果\n",
    "median_blur_result_5_5 = cv.medianBlur(src_img, 5)\n",
    "my_plot(2, 2, 4, median_blur_result_5_5, \"Median Blur(5*5)\")\n",
    "\n",
    "plt.show()\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果:\n",
    "![原图与中值滤波结果的比较](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/original_median_blur.png)\n",
    "#### 结果分析：\n",
    "+ 3*3的中值滤波结果与原图相比，平滑效果很明显，特别是脸部，微小的皱纹不见了，又极大程度上保留了原图的信息；\n",
    "+ 5*5的中值滤波结果与原图相比，平滑效果达到了，但变得较为模糊；"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (4) 在两种窗口大小的情况下，平均滤波、高斯滤波和中值滤波的效果比较"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "def my_plot(rows, cols, seq, img, title):\n",
    "    plt.subplot(rows, cols, seq)\n",
    "    plt.title(title)\n",
    "    plt.imshow(reverse(img))\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "\n",
    "    \n",
    "src_img = cv.imread(r\"image/lena.jpg\")\n",
    "plt.figure(figsize=(2.5 * src_img.shape[0] * 0.01, 2.5 * src_img.shape[1] * 0.01))\n",
    "my_plot(3, 3, 1, src_img, \"Original\")\n",
    "\n",
    "# 3*3窗口平均滤波结果\n",
    "averaging_blur_result_3_3 = cv.blur(src_img, (3, 3))\n",
    "my_plot(3, 3, 4, averaging_blur_result_3_3, \"Averaging Blur(3*3)\")\n",
    "\n",
    "# 3*3窗口高斯滤波结果\n",
    "gaussian_blur_result_3_3 = cv.GaussianBlur(src_img, (3, 3), 0)\n",
    "my_plot(3, 3, 5, gaussian_blur_result_3_3, \"Gaussian Blur(3*3)\")\n",
    "\n",
    "# 3*3窗口中值滤波结果\n",
    "median_blur_result_3_3 = cv.medianBlur(src_img, 3)\n",
    "my_plot(3, 3, 6, median_blur_result_3_3, \"Median Blur(3*3)\")\n",
    "\n",
    "# 5*5窗口平均滤波结果\n",
    "averaging_blur_result_5_5 = cv.blur(src_img, (5, 5))\n",
    "my_plot(3, 3, 7, averaging_blur_result_5_5, \"Averaging Blur(5*5)\")\n",
    "\n",
    "# 5*5窗口高斯滤波结果\n",
    "gaussian_blur_result_5_5 = cv.GaussianBlur(src_img, (5, 5), 0)\n",
    "my_plot(3, 3, 8, gaussian_blur_result_5_5, \"Gaussian Blur(5*5)\")\n",
    "\n",
    "# 5*5窗口中值滤波结果\n",
    "median_blur_result_5_5 = cv.medianBlur(src_img, 5)\n",
    "my_plot(3, 3, 9, median_blur_result_5_5, \"Median Blur(5*5)\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果:\n",
    "![原图与3*3的平均滤波、高斯滤波和中值滤波结果比较](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/original_averaging_gaussian_median.png)\n",
    "#### 结果分析:\n",
    "+ 3\\*3窗口的平均滤波、高斯滤波、中值滤波，都对原图有一定程度的平滑效果，且都保留了原图的信息;\n",
    "+ 3\\*3窗口的高斯滤波和中值滤波的效果相当，它们稍好于平均滤波;\n",
    "+ 当窗口尺寸变为5\\*5时，可以看到平均滤波的效果像毛玻璃的感觉，已经有点模糊了；但高斯和中值的结果保持的不错。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (5) 用一张椒盐噪声较多的图像作为原图，来测试平均滤波、高斯滤波和中值滤波的效果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "def my_plot(rows, cols, seq, img, title):\n",
    "    plt.subplot(rows, cols, seq)\n",
    "    plt.title(title)\n",
    "    plt.imshow(reverse(img))\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "\n",
    "\n",
    "# 显示原图\n",
    "src_img = cv.imread(r\"image/pic2.png\")\n",
    "plt.figure(2, (3 * src_img.shape[0] * 0.01, 4 * src_img.shape[1] * 0.01))\n",
    "my_plot(4, 3, 1, src_img, \"Original\")\n",
    "\n",
    "# 显示灰度图\n",
    "src_img_gray = cv.imread(r\"image/pic2.png\", 0)\n",
    "plt.subplot(4, 3, 2), plt.title(\"Original gray\"), \\\n",
    "    plt.imshow(src_img_gray, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "# 显示灰度直方图\n",
    "hist_cv = cv.calcHist([src_img_gray], [0], None, [256], [0, 256])\n",
    "plt.subplot(4, 3, 3), plt.title(\"Gray Histogram\"), plt.plot(hist_cv)\n",
    "\n",
    "# 3*3窗口平均滤波\n",
    "averaging_blur_result_3_3 = cv.blur(src_img, (3, 3))\n",
    "my_plot(4, 3, 4, averaging_blur_result_3_3, \"Averaging Blur(3*3)\")\n",
    "\n",
    "# 3*3窗口高斯滤波\n",
    "gaussian_blur_result_3_3 = cv.GaussianBlur(src_img, (3, 3), 0)\n",
    "my_plot(4, 3, 5, gaussian_blur_result_3_3, \"Gaussian Blur(3*3)\")\n",
    "\n",
    "# 3*3窗口中值滤波\n",
    "median_blur_result_3_3 = cv.medianBlur(src_img, 3)\n",
    "my_plot(4, 3, 6, median_blur_result_3_3, \"Median Blur(3*3)\")\n",
    "\n",
    "# 7*7窗口平均滤波\n",
    "averaging_blur_result_7_7 = cv.blur(src_img, (7, 7))\n",
    "my_plot(4, 3, 7, averaging_blur_result_7_7, \"Averaging Blur(7*7)\")\n",
    "\n",
    "# 7*7窗口高斯滤波\n",
    "gaussian_blur_result_7_7 = cv.GaussianBlur(src_img, (7, 7), 0)\n",
    "my_plot(4, 3, 8, gaussian_blur_result_7_7, \"Gaussian Blur(7*7)\")\n",
    "\n",
    "# 7*7窗口中值滤波\n",
    "median_blur_result_7_7 = cv.medianBlur(src_img, 7)\n",
    "my_plot(4, 3, 9, median_blur_result_7_7, \"Median Blur(7*7)\")\n",
    "\n",
    "# 7*7窗口中值滤波后的灰度图和灰度直方图\n",
    "median_blur_result_7_7_gray = cv.cvtColor(averaging_blur_result_7_7, \n",
    "                                          cv.COLOR_BGR2GRAY)\n",
    "plt.subplot(4, 3, 10), \\\n",
    "    plt.title(\"Median Blur(7*7) Gray\"), \\\n",
    "    plt.imshow(median_blur_result_7_7_gray, \"gray\"),\\\n",
    "    plt.xticks([]), plt.yticks([])\n",
    "median_blur_result_7_7_gray_histo = cv.calcHist([median_blur_result_7_7_gray], \n",
    "                                                [0], None, [256], [0, 256])\n",
    "plt.subplot(4, 3, 11), plt.title(\"Median Blur(7*7) Gray Histogram\"), \\\n",
    "    plt.plot(median_blur_result_7_7_gray_histo)\n",
    "\n",
    "plt.show()\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果：\n",
    "![平均-高斯-中值滤波对椒盐的处理效果对比](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/original_three_blur.png)\n",
    "#### 结果分析：\n",
    "+ 通过查看原图的灰度直方图，可以看出，有一些灰度的小尖峰，根据原图特点判断，应该是属于椒盐噪声导致；\n",
    "+ 通过看第二行和第三行的滤波效果图，对应一张含有较多椒盐噪声的原图，在3\\*3的窗口下和7\\*7的窗口下，中值滤波的去噪效果最好，平均和高斯滤波的去噪效果差不多。本质的原因，还是因为中值滤波的目标图像的像素值，是从原图中获取的，而高斯和平均滤波，是采用了外部确定的加权系数；因此对于这种椒盐噪声较多的图像，是中值滤波去噪更为有效。\n",
    "+ 第四行，进一步的，得到中值滤波的结果灰度图和灰度直方图，通过灰度直方图可以看出，灰度值有两个峰值点，一个是在100附近，一个是在200附近，分别对应三个小黑块，和其余偏白的部分，没有再出现多余的小尖峰，可见去噪效果不错。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (6) 分析双边滤波的效果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "def my_plot(rows, cols, seq, img, title):\n",
    "    plt.subplot(rows, cols, seq)\n",
    "    plt.title(title)\n",
    "    plt.imshow(reverse(img))\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "\n",
    "\n",
    "src_img = cv.imread(r\"image/bilateral.jpg\")\n",
    "\n",
    "plt.figure(figsize=(2.5 * src_img.shape[0] * 0.01, 0.8 * src_img.shape[1] * 0.01))\n",
    "\n",
    "my_plot(1, 2, 1, src_img, \"Original\")\n",
    "\n",
    "blurred_img = cv.bilateralFilter(src_img, 9, 75, 75)\n",
    "my_plot(1, 2, 2, blurred_img, \"BilateralFiltered Image\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果\n",
    "![双边滤波结果](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/original_bilateral_blur.png)\n",
    "#### 结果分析\n",
    "+ 可以看到，图像的纹理部分得到了平滑，但所有边缘都保留下来了。\n",
    "+ 进一步研究双边滤波的原理，可以知道，它的处理过程中，使用了两个高斯函数，一个用来确定卷积核的元素数值，一个用来表达差分值的分布；前者的作用是对指定区域进行高斯模糊，后者的作用是根据差分值阈值的判断，差分值超过阈值的像素，被认为是边缘部分，不参与高斯模糊的处理。\n",
    "+ 这种机制的结果是，既保留了边缘，又平滑了非边缘部分，以便图像的后续处理。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. 以Lena为原始图像，通过OpenCV使用Sobel及Canny算子检测，比较边缘检测结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "def my_plot(rows, cols, seq, img, title):\n",
    "    plt.subplot(rows, cols, seq)\n",
    "    plt.title(title)\n",
    "    plt.imshow(reverse(img))\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "\n",
    "\n",
    "# 画出原图\n",
    "src_img = cv.imread(r\"image/lena.jpg\")\n",
    "plt.figure(figsize=(3.5 * src_img.shape[0] * 0.01, 3.5 * src_img.shape[1] * 0.01))\n",
    "my_plot(3, 3, 1, src_img, \"Original\")\n",
    "\n",
    "# 画出Sobel变换效果图\n",
    "sobel_x = cv.Sobel(src_img, cv.CV_16S, 1, 0)\n",
    "sobel_x_8u = cv.convertScaleAbs(sobel_x)\n",
    "my_plot(3, 3, 4, sobel_x_8u, \"Sobel X\")\n",
    "\n",
    "sobel_y = cv.Sobel(src_img, cv.CV_16S, 0, 1)\n",
    "sobel_y_8u = cv.convertScaleAbs(sobel_y)\n",
    "my_plot(3, 3, 5, sobel_y_8u, \"Sobel Y\")\n",
    "\n",
    "sobel = cv.addWeighted(sobel_x_8u, 0.5, sobel_y_8u, 0.5, 0)\n",
    "my_plot(3, 3, 6, sobel, \"Sobel\")\n",
    "\n",
    "# 画出Canny变换效果图\n",
    "src_img_gray = cv.imread(r\"image/lena.jpg\", 0)\n",
    "src_img_blur = cv.GaussianBlur(src_img_gray, (3, 3), 0)\n",
    "canny = cv.Canny(src_img_blur, 50, 160)\n",
    "plt.subplot(3, 3, 7), plt.imshow(canny, \"gray\"), \\\n",
    "plt.title(\"Canny(minVal=15,maxVal=30)\"), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果\n",
    "![Sobel和Canny边缘检测效果比较](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/sobel_canny.png)\n",
    "#### 结果分析\n",
    "+ Sobel变换检测出来的边缘，存在两个问题\n",
    "\t+ 误检：边缘线较粗，不够精细；\n",
    "\t+ 漏检：边缘线看起来断断续续，轮廓不是很完整；\n",
    "**原因分析**：Sobel的思想，是对每个像素求卷积，通过卷积核的设计，卷积的过程的实质是在求x,y两个方向的图像差分值，来近似模拟梯度值；然后对新生成的像素灰度值做阈值处理，处理后的目标图像，就只剩下了边缘部分；实际运算过程中，为了计算简便，用x,y方向的差分值的绝对值之和来近似当做梯度值；正因为计算方向的单一性，遇到灰度值分布稍微复杂，边缘处的像素差分值并不都满足统一的阈值的情况，就容易发生误检(边缘较粗)或者漏检(边缘线断续)的情况，所以，它的适用对象是边缘灰度值变化较为统一的图像。\n",
    "+ Canny算子检测的结果\n",
    "\t+ 好的信噪比，即将非边缘判定为边缘的概率要低：这个基本达到，依赖于高斯模糊的预处理，以及将maxVal设置大一些;\n",
    "\t+ 高的定位性能，即检测出的边缘点要尽可能在实际边缘的中心：这个也达到了，但是仍然受原图亮度的影响较大，导致有些较为暗的边缘部分，没有被当做中心检测出来，比如帽檐；\n",
    "\t+ 对单一边缘仅有唯一的响应，即单个边缘产生多个响应的概率要低：这个基本满足；\n",
    "\n",
    "+ Canny和Sobel的比较\n",
    "\t+ Canny算子只能处理灰度图，但彩色图像转灰度图，通常有10%左右的边缘损失，进而会影响检测的准确性。\n",
    "\t+ Canny算子需要根据图像情况手动调整阈值，不过针对自适应调整阈值的方向，已经有了一些改进，比如小波变换、信息熵、BP神经网络、灰度矩等方法。\n",
    "\t+ Canny算子的检测结果连通性较好，断裂的边缘相比Sobel算子的结果大大减少，这个取决于minVal的设置，设置过大会导致断裂的情况增多；\n",
    "\t+ Sobel基于一阶导数，所以检测出的边缘稍微粗一些，而Canny算子虽然也是基于一阶导数，但因为增加了非极大值抑制的步骤，所以边缘精细了很多。\n",
    "\t+ 处理速度方面，Sobel因为简单，所以快一些；"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. 在OpenCV安装目录下找到课程对应演示图片(安装目录\\sources\\samples\\data)，首先计算灰度直方图，进一步使用大津算法进行分割，并比较分析分割结果。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import copy\n",
    "\n",
    "src_img = cv.imread(r\"image/rice.png\", 0)\n",
    "\n",
    "plt.figure(figsize=(3.5 * src_img.shape[0] * 0.01, 5.5 * src_img.shape[1] * 0.01))\n",
    "\n",
    "\n",
    "def deal_original(src_img):\n",
    "    \"\"\"\n",
    "    处理原图(灰度直方图+大津处理)\n",
    "    :param src_img: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    # 画出原图的灰度图\n",
    "    plt.subplot(5, 3, 1), plt.title(\"Original Gray\"), plt.imshow(src_img, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "    # 原图的灰度直方图\n",
    "    src_gray_histo = cv.calcHist([src_img], [0], None, [256], [0, 256])\n",
    "    plt.subplot(5, 3, 2), plt.title(\"Original Gray Histogram\"), \\\n",
    "        plt.plot(src_gray_histo), \\\n",
    "        plt.xlabel(\"Gray Scale\"), plt.ylabel(\"Pixel Num\"), plt.xlim([0, 256])\n",
    "    # 原图经过大津处理的结果\n",
    "    threshold, src_otsu = cv.threshold(src_img, 0, 255, cv.THRESH_OTSU)\n",
    "    plt.subplot(5, 3, 3), plt.title(\"Original OTSU[threshold:{}]\".format(threshold)), \\\n",
    "        plt.imshow(src_otsu, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def deal_gaussian(src_img):\n",
    "    \"\"\"\n",
    "    先对原图进行高斯模糊预处理降噪,再画出灰度直方图和大津处理效果图\n",
    "    :param src_img: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    gaussian_blur = cv.GaussianBlur(src_img, (3, 3), 0)\n",
    "    plt.subplot(5, 3, 4), plt.title(\"Gaussian Blur(3*3)\"), plt.imshow(gaussian_blur, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "    \n",
    "    gaussian_gray_histo = cv.calcHist([gaussian_blur], [0], None, [256], [0, 256])\n",
    "    plt.subplot(5, 3, 5), plt.title(\"Gaussian Gray Histogram\"), \\\n",
    "        plt.plot(gaussian_gray_histo),\\\n",
    "        plt.xlabel(\"Gray Scale\"), plt.ylabel(\"Pixel Num\"), plt.xlim([0, 256])\n",
    "    \n",
    "    threshold, gaussian_otsu = cv.threshold(gaussian_blur, 0, 255, cv.THRESH_OTSU)\n",
    "    plt.subplot(5, 3, 6), plt.title(\"Gaussian OTSU[threshold:{}]\".format(threshold)),\\\n",
    "        plt.imshow(gaussian_otsu, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def deal_morphological(src_img):\n",
    "    \"\"\"\n",
    "    先对原图进行形态学预处理降噪,再画出灰度直方图和大津处理效果图\n",
    "    :param src_img: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    kernel = np.ones((3, 3), np.uint8)\n",
    "    erosion = cv.erode(src_img, kernel)\n",
    "    dilation = cv.dilate(erosion, kernel, iterations=1)\n",
    "    plt.subplot(5, 3, 7), plt.title(\"Erosion then Dilation\"), plt.imshow(dilation, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "    \n",
    "    morphological_gray_histo = cv.calcHist([dilation], [0], None, [256], [0, 256])\n",
    "    plt.subplot(5, 3, 8), plt.title(\"Morphological Gray Histogram\"), \\\n",
    "        plt.plot(morphological_gray_histo),\\\n",
    "        plt.xlabel(\"Gray Scale\"), plt.ylabel(\"Pixel Num\"), plt.xlim([0, 256])\n",
    "\n",
    "    threshold, morphological_otsu = cv.threshold(dilation, 0, 255, cv.THRESH_OTSU)\n",
    "    plt.subplot(5, 3, 9), plt.title(\"Morphological OTSU[threshold:{}]\".format(threshold)), \\\n",
    "        plt.imshow(morphological_otsu, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def deal_gaussian_adaptive_threshold(src_img):\n",
    "    \"\"\"\n",
    "    先对原图进行高斯模糊预处理降噪,再画出灰度直方图和局部阈值分割处理效果图\n",
    "    :param src_img: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    gaussian_blur = cv.GaussianBlur(src_img, (3, 3), 0)\n",
    "    plt.subplot(5, 3, 10), plt.title(\"Gaussian Blur(3*3)\"), plt.imshow(gaussian_blur, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "    \n",
    "    gaussian_gray_histo = cv.calcHist([gaussian_blur], [0], None, [256], [0, 256])\n",
    "    plt.subplot(5, 3, 11), plt.title(\"Gaussian Gray Histogram\"), \\\n",
    "        plt.plot(gaussian_gray_histo),\\\n",
    "        plt.xlabel(\"Gray Scale\"), plt.ylabel(\"Pixel Num\"), plt.xlim([0, 256])\n",
    "    \n",
    "    gaussian_adaptive_threold = cv.adaptiveThreshold(dilation, 255, cv.ADAPTIVE_THRESH_MEAN_C, cv.THRESH_BINARY, 51, 1)\n",
    "    plt.subplot(5, 3, 12), plt.title(\"Gaussian Adaptive Threshold\"), \\\n",
    "        plt.imshow(gaussian_adaptive_threold, \"gray\"),\\\n",
    "        plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def deal_morphological_adaptive_threshold(src_img):\n",
    "    \"\"\"\n",
    "    先对原图进行形态学预处理降噪,再画出灰度直方图和局部阈值分割处理效果图\n",
    "    :param src_img: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    kernel = np.ones((3, 3), np.uint8)\n",
    "    erosion = cv.erode(src_img, kernel)\n",
    "    dilation = cv.dilate(erosion, kernel, iterations=1)\n",
    "    plt.subplot(5, 3, 13), plt.title(\"Erosion then Dilation\"), plt.imshow(dilation, \"gray\"), plt.xticks([]), plt.yticks([])\n",
    "    \n",
    "    morphological_gray_histo = cv.calcHist([dilation], [0], None, [256], [0, 256])\n",
    "    plt.subplot(5, 3, 14), plt.title(\"Morphological Gray Histogram\"), \\\n",
    "        plt.plot(morphological_gray_histo),\\\n",
    "        plt.xlabel(\"Gray Scale\"), plt.ylabel(\"Pixel Num\"), plt.xlim([0, 256])\n",
    "    \n",
    "    morphological_adaptive_threold = cv.adaptiveThreshold(dilation, 255, cv.ADAPTIVE_THRESH_MEAN_C, cv.THRESH_BINARY, 51, 1)\n",
    "    plt.subplot(5, 3, 15), plt.title(\"Morphological Adaptive Threshold\"), \\\n",
    "        plt.imshow(morphological_adaptive_threold, \"gray\"),\\\n",
    "        plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "src_img_copy_1 = copy.deepcopy(src_img)\n",
    "deal_original(src_img_copy_1)\n",
    "\n",
    "src_img_copy_2 = copy.deepcopy(src_img)\n",
    "deal_gaussian(src_img_copy_2)\n",
    "\n",
    "src_img_copy_3 = copy.deepcopy(src_img)\n",
    "deal_morphological(src_img_copy_3)\n",
    "\n",
    "src_img_copy_4 = copy.deepcopy(src_img)\n",
    "deal_gaussian_adaptive_threshold(src_img_copy_4)\n",
    "\n",
    "src_img_copy_5 = copy.deepcopy(src_img)\n",
    "deal_morphological_adaptive_threshold(src_img_copy_5)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果\n",
    "![OTSU结果](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/otsu.png)\n",
    "#### 结果分析\n",
    "+ 由原图的灰度直方图可以看出，灰度主要分布三个尖峰值附近,\n",
    "\t+ 其中灰度值最大的那个尖峰，对应白色米粒；\n",
    "\t+ 中间的尖峰，对应原图上面大部分的灰灰的背景；\n",
    "\t+ 灰度值最小的尖峰，对应图像下方的偏暗的背景；\n",
    "+ 对原图直接进行OTSU阈值化，得到的结果含有很多噪声，得到的OTSU阈值也偏高；\n",
    "+ 对原图先进行高斯模糊的预处理，再进行OTSU阈值化，对噪声的去除效果较好，但残存少部分噪声点，得到的OTSU阈值居中；\n",
    "+ 对原图进行先腐蚀再膨胀的预处理，再进行OTSU阈值化，对噪声的去除效果是最好的，但也误伤了部分米粒，得到的OTSU阈值偏低。\n",
    "\n",
    "#### 进一步处理\n",
    "+ 鉴于原图下方偏暗，所以采用全局统一的阈值，可能难免会导致部分米粒阈值化后变残，所以需要尽可能用到局部特性，所以选择局部阈值化分割方法来进行处理\n",
    "+ 通过对比看出，不管选择形态学还是高斯模糊先去噪，再进行局部阈值分割，分割后的米粒跟原图中的大小基本保持一致，而且噪声去除效果也很好。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. 使用米粒图像，分割得到各米粒，首先计算各区域(米粒)的面积、长度等信息，进一步计算面积、长度的均值及方差，分析落在3sigma范围内米粒的数量。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import copy\n",
    "\n",
    "src_img = cv.imread(r\"image/rice.png\")\n",
    "\n",
    "src_img_copy = copy.deepcopy(src_img)\n",
    "\n",
    "src_img = cv.cvtColor(src_img, cv.COLOR_BGR2GRAY)\n",
    "\n",
    "plt.figure(figsize=(3.5 * src_img.shape[0] * 0.01, 2.5 * src_img.shape[1] * 0.01))\n",
    "\n",
    "# 画出原图\n",
    "plt.subplot(2, 3, 1), plt.title(\"Original Gray\"), plt.imshow(src_img, \"gray\"), \\\n",
    "    plt.xticks([]), plt.yticks([])\n",
    "\n",
    "# 进行形态学处理\n",
    "kernel = np.ones((3, 3), np.uint8)\n",
    "result = cv.erode(src_img, kernel, iterations=2)\n",
    "result = cv.dilate(result, kernel, iterations=2)\n",
    "plt.subplot(2, 3, 2), plt.title(\"Morphological\"), plt.imshow(result, \"gray\"), \\\n",
    "    plt.xticks([]), plt.yticks([])\n",
    "\n",
    "# 进行局部阈值分割\n",
    "result = cv.adaptiveThreshold(result, 255, cv.ADAPTIVE_THRESH_MEAN_C, cv.THRESH_BINARY, 61, 1)\n",
    "plt.subplot(2, 3, 3), plt.title(\"Adaptive Threshold\"), plt.imshow(result, \"gray\"), \\\n",
    "    plt.xticks([]), plt.yticks([])\n",
    "\n",
    "# 进行轮廓提取\n",
    "image, contours, hierarchy = cv.findContours(result, cv.RETR_TREE, cv.CHAIN_APPROX_SIMPLE)\n",
    "pre_contour_copy = copy.deepcopy(src_img_copy)\n",
    "contour_image = cv.drawContours(pre_contour_copy, contours, -1, (255, 0, 0), 1)\n",
    "plt.subplot(2, 3, 4), plt.title(\"Contours\"), plt.imshow(contour_image), \\\n",
    "    plt.xticks([]), plt.yticks([])\n",
    "\n",
    "# 计算各个米粒的面积、长度(最小包围矩形的长)\n",
    "areas = []\n",
    "lens = []\n",
    "for i in range(1, len(contours) + 1):\n",
    "    contour = contours[i - 1]\n",
    "    area = cv.contourArea(contour)\n",
    "    \n",
    "    if area < 10:\n",
    "        continue\n",
    "    \n",
    "    areas.append(int(area))\n",
    "    \n",
    "    minRect = cv.minAreaRect(contour)   # ((center_x, center_y), (width, height), rotation_angle)\n",
    "    box = cv.boxPoints(minRect)         # [[x1, y1], [x2, y2], [x3, y3], [x4, y4]]\n",
    "    box = np.int0(box)                  # [[x1, y1], [x2, y2], [x3, y3], [x4, y4]]\n",
    "    \n",
    "    cv.drawContours(src_img_copy, [box], 0, (0, 255, 0), 1)\n",
    "    cv.putText(src_img_copy, str(len(areas)), tuple(box[1]), cv.FONT_HERSHEY_PLAIN, 0.5, (255, 0, 0))\n",
    "    \n",
    "    lens.append(int(max(minRect[1])))\n",
    "\n",
    "plt.subplot(2, 3, 5), plt.title(\"Contours Rect\"), plt.imshow(src_img_copy, \"gray\"), \\\n",
    "    plt.xticks([]), plt.yticks([])\n",
    "\n",
    "count = len(areas)\n",
    "print(\"检测出的米粒数:{}\".format(count))\n",
    "print(\"\\n米粒的面积:{}\".format(areas))\n",
    "print(\"\\n米粒的长度:{}\".format(lens))\n",
    "\n",
    "# 计算面积的均值和标准差\n",
    "area_mean = int(np.mean(areas))\n",
    "area_sigma = int(np.std(areas, ddof=1))\n",
    "area_mean_min, area_mean_max = area_mean - 3 * area_sigma, area_mean + 3 * area_sigma\n",
    "area_best_num = len(list(filter(lambda x: area_mean_min <= x <= area_mean_max, areas)))\n",
    "print(\"\\n米粒的面积均值: {}, 标准差: {}, 面积落在[{} - 3 * {}, {} + 3 * {}]区域的米粒数量: {}\".format(area_mean, \n",
    "                                                            area_sigma, area_mean, area_sigma, area_mean, area_sigma,\n",
    "                                                            area_best_num))\n",
    "\n",
    "# 计算长度的均值和标准差\n",
    "len_mean = int(np.mean(lens))\n",
    "len_sigma = int(np.std(lens, ddof=1))\n",
    "len_mean_min, len_mean_max = len_mean - 3 * len_sigma, len_mean + 3 * len_sigma\n",
    "len_best_num = len(list(filter(lambda x: len_mean_min <= x <= len_mean_max, lens)))\n",
    "print(\"\\n米粒的长度均值: {}, 标准差: {}, \"\n",
    "      \"长度落在[{} - 3 * {}, {} + 3 * {}]\"\n",
    "      \"区域的米粒数量: {}\".format(len_mean,\n",
    "                           len_sigma, \n",
    "                           len_mean, \n",
    "                           len_sigma, \n",
    "                           len_mean, \n",
    "                           len_sigma,\n",
    "                           len_best_num))\n",
    "\n",
    "plt.show()\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 运行结果\n",
    "![米粒检测结果图](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/rice_detect.png)\n",
    "\n",
    "检测出的米粒数:95\n",
    "\n",
    "米粒的面积:[132, 155, 204, 180, 256, 231, 138, 122, 93, 197, 209, 131, 179, 210, 222, 189, 186, 172, 203, 213, 188, 68, 206, 197, 165, 215, 189, 265, 109, 182, 142, 202, 183, 171, 201, 129, 204, 160, 212, 212, 150, 176, 164, 152, 214, 197, 145, 143, 165, 139, 147, 204, 198, 214, 181, 196, 34, 91, 352, 376, 188, 172, 173, 160, 98, 140, 189, 167, 117, 158, 128, 175, 176, 182, 200, 153, 141, 168, 114, 162, 369, 189, 169, 162, 126, 179, 98, 66, 35, 154, 101, 12, 170, 105, 48]\n",
    "\n",
    "米粒的长度:[20, 21, 26, 28, 28, 29, 21, 17, 19, 27, 28, 21, 23, 26, 29, 25, 28, 23, 28, 28, 24, 11, 28, 27, 25, 29, 28, 37, 18, 26, 24, 26, 28, 24, 29, 20, 28, 21, 27, 28, 23, 24, 22, 24, 27, 26, 22, 24, 25, 24, 25, 27, 27, 30, 27, 26, 16, 12, 36, 38, 26, 27, 25, 25, 16, 21, 28, 27, 15, 24, 21, 25, 26, 27, 27, 23, 21, 24, 21, 24, 42, 26, 24, 24, 19, 26, 16, 12, 8, 23, 18, 6, 24, 17, 8]\n",
    "\n",
    "米粒的面积均值: 168, 标准差: 58, 面积落在[168 - 3 * 58, 168 + 3 * 58]区域的米粒数量: 92\n",
    "\n",
    "米粒的长度均值: 23, 标准差: 5, 长度落在[23 - 3 * 5, 23 + 3 * 5]区域的米粒数量: 93\n",
    "\n",
    "#### 结果分析\n",
    "+ 先通过形态学处理(先2次腐蚀, 再2次膨胀)，可以看到第2张图噪声已被滤除，且米粒边缘也更加简洁；\n",
    "+ 再通过区域阈值分割，将灰度图转化为二值图像；\n",
    "+ 接着通过轮廓提取，画出了轮廓图，可以看到包含一些较小的米粒在内了；\n",
    "+ 然后根据上一步提取的轮廓点，在原图的复制品上，画出每个米粒的最小包围矩形和序号，并记录每个米粒轮廓的面积，和最小包围矩形的长度；\n",
    "可以看到有3处较为接近的2个米粒连在了一起；\n",
    "+ 计算米粒面积的均值和方差，以及落在3sigma范围内的米粒数：共计95个米粒，面积均值168，标准差58，面积落在3sigma范围内的米粒数为92，占比96.8%;\n",
    "+ 计算米粒长度的均值和方差，以及落在3sigma范围内的米粒数：攻击95个米粒，长度均值为23，标准差5，长度落在3sigma范围内的米粒数为93，占比97.9%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. 使用棋盘格及自选风景图像，分别使用SIFT、FAST及ORB算子检测角点，并比较分析检测结果。 \n",
    "### (可选)使用Harris角点检测算子检测棋盘格，并与上述结果比较。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding:utf8 -*-\n",
    "import cv2 as cv\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import copy\n",
    "\n",
    "\n",
    "def reverse(img):\n",
    "    \"\"\"\n",
    "    将BGR转化为RGB\n",
    "    :param img: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    b, g, r = cv.split(img)\n",
    "    return cv.merge([r, g, b])\n",
    "\n",
    "\n",
    "# 读取棋盘格图像\n",
    "# src_img = cv.imread(r\"image/chessboard.png\")\n",
    "\n",
    "# 读取风景图像\n",
    "src_img = cv.imread(r\"image/scene.jpg\")\n",
    "\n",
    "# 设置画布\n",
    "plt.figure(figsize=(2.1 * src_img.shape[0] * 0.01, 1.1 * src_img.shape[1] * 0.01))\n",
    "\n",
    "# 画出原图\n",
    "plt.subplot(1, 2, 1), plt.title(\"Original\"), plt.imshow(reverse(src_img)),\\\n",
    "    plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def cal_histo(src_img_copy):\n",
    "    \"\"\"\n",
    "    画出灰度直方图\n",
    "    :param src_img_copy: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    src_img_copy_gray = cv.cvtColor(src_img_copy, cv.COLOR_BGR2GRAY)\n",
    "    histogram = cv.calcHist([src_img_copy_gray], [0], None, [256], [0, 256])\n",
    "    plt.subplot(1, 2, 2), plt.title(\"Histogram\"), plt.plot(histogram),\\\n",
    "        plt.xlabel(\"Gray Scale\"), plt.ylabel(\"Pixel Num\")\n",
    "\n",
    "\n",
    "def harris_detect(src_img_copy):\n",
    "    \"\"\"\n",
    "    执行Harris角点检测\n",
    "    :param src_img_copy: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    src_img_gray = np.float32(cv.cvtColor(src_img_copy, cv.COLOR_BGR2GRAY))\n",
    "    dst = cv.cornerHarris(src_img_gray, 2, 3, 0.04)\n",
    "    dst = cv.dilate(dst, None)\n",
    "    dst = np.int0(dst)\n",
    "    \n",
    "    dst_zero_threshold = 0.01 * dst.max()\n",
    "    for i in range(0, dst.shape[0]):\n",
    "        for j in range(0, dst.shape[1]):\n",
    "            if dst[i][j] > dst_zero_threshold:\n",
    "                cv.circle(src_img_copy, (j, i), 1, (0, 0, 255), -1)\n",
    "    plt.subplot(1, 2, 2), plt.title(\"Harris\"), plt.imshow(reverse(src_img_copy)),\\\n",
    "        plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def sift_detect(src_img_copy):\n",
    "    \"\"\"\n",
    "    SIFT检测角点\n",
    "    :param src_img_copy: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    src_img_copy_gray = cv.cvtColor(src_img_copy, cv.COLOR_BGR2GRAY)\n",
    "    # pixel_max, pixel_min = np.amax(src_img_copy_gray), np.amin(src_img_copy_gray)\n",
    "    n_octave_layers = 3\n",
    "    # contrast_threshold = 0.01 * n_octave_layers * (pixel_max - pixel_min) / 255.0\n",
    "    contrast_threshold = 0.04\n",
    "    sift = cv.xfeatures2d.SIFT_create(0, n_octave_layers, contrast_threshold,\n",
    "                                      10, 1.6)\n",
    "    # kps = sift.detect(src_img_copy_gray)\n",
    "    # features = sift.compute(src_img_copy_gray, kps)\n",
    "    \n",
    "    kps, features = sift.detectAndCompute(src_img_copy_gray, None)\n",
    "    \n",
    "    for kp in kps:\n",
    "        x, y = int(kp.pt[0]), int(kp.pt[1])\n",
    "        color = (0, 0, 255)\n",
    "        radius = 3\n",
    "        cv.circle(src_img_copy, (x, y), radius, color, -1)\n",
    "    plt.subplot(1, 2, 2), plt.title(\"SIFT\"), plt.imshow(reverse(src_img_copy)), \\\n",
    "        plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def fast_detect(src_img_copy):\n",
    "    \"\"\"\n",
    "    FAST检测角点\n",
    "    :param src_img_copy: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    src_img_copy_gray = cv.cvtColor(src_img_copy, cv.COLOR_BGR2GRAY)\n",
    "    fast = cv.FastFeatureDetector_create(nonmaxSuppression=True, \n",
    "                                         type=cv.FAST_FEATURE_DETECTOR_TYPE_7_12)\n",
    "    kps = fast.detect(src_img_copy_gray, None)\n",
    "    for kp in kps:\n",
    "        x, y = int(kp.pt[0]), int(kp.pt[1])\n",
    "        color = (0, 0, 255)\n",
    "        radius = 1\n",
    "        cv.circle(src_img_copy, (x, y), radius, color, -1)\n",
    "    \n",
    "    print(\"type:{},nonMaxSuppression:{},threshold={}\".format(\n",
    "        fast.getType(),\n",
    "        fast.getNonmaxSuppression(),\n",
    "        fast.getThreshold()\n",
    "    ))\n",
    "    \n",
    "    plt.subplot(1, 2, 2), plt.title(\"FAST\"), plt.imshow(reverse(src_img_copy)),\\\n",
    "        plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "def orb_detect(src_img_copy):\n",
    "    \"\"\"\n",
    "    ORB检测角点\n",
    "    :param src_img_copy: \n",
    "    :return: \n",
    "    \"\"\"\n",
    "    src_img_copy_gray = cv.cvtColor(src_img_copy, cv.COLOR_BGR2GRAY)\n",
    "    orb = cv.ORB_create()\n",
    "    kps = orb.detect(src_img_copy_gray, None)\n",
    "    for kp in kps:\n",
    "        x, y = int(kp.pt[0]), int(kp.pt[1])\n",
    "        radius = 1\n",
    "        color = (0, 0, 255)\n",
    "        cv.circle(src_img_copy, (x, y), radius, color, -1)\n",
    "    \n",
    "    plt.subplot(1, 2, 2), plt.title(\"ORB\"), plt.imshow(reverse(src_img_copy)),\\\n",
    "        plt.xticks([]), plt.yticks([])\n",
    "\n",
    "\n",
    "# cal_histo(copy.deepcopy(src_img))\n",
    "# harris_detect(copy.deepcopy(src_img))\n",
    "# sift_detect(copy.deepcopy(src_img))\n",
    "# fast_detect(copy.deepcopy(src_img))\n",
    "# orb_detect(copy.deepcopy(src_img))\n",
    "\n",
    "\n",
    "plt.show()\n",
    "plt.waitforbuttonpress()\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 对于棋盘格的运行结果\n",
    "+ Harris角点检测\n",
    "![Harris角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/harris.png)\n",
    "\n",
    "+ SIFT角点检测\n",
    "全部默认参数情形：\n",
    "![SIFT角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/sift.png)\n",
    "nfeatures=100\n",
    "![SIFT角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/sift_nfeatures_100.png)\n",
    "nfeatures=124\n",
    "![SIFT角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/sift_nfeatures_124.png)\n",
    "nfeatures=500\n",
    "![SIFT角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/sift_nfeatures_500.png)\n",
    "nfeatures=1000\n",
    "![SIFT角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/sift_nfeatures_1000.png)\n",
    "\n",
    "+ FAST角点检测\n",
    "![FAST角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/fast.png)\n",
    "\n",
    "+ ORB角点检测\n",
    "![ORB角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/orb.png)\n",
    "\n",
    "\n",
    "#### 对于棋盘格的检测结果分析\n",
    "+ 灰度直方图显示，图中只含有灰度值为0和255的两种像素，分别对应黑白格；\n",
    "+ Harris角点检测特别适合这种角点处对比特别明显的图，角点检测结果十分精准，运行速度很快；\n",
    "+ SIFT角点检测的结果，因为要照顾到多尺度的情形，所以检测的特征点往图右略偏，左边一些角点未被检测出来；运行速度比较慢；\n",
    "实验过程中尝试了SIFT参数的一些调整，其中\n",
    "\t+ SIFT_create(nfeatures=None, nOctaveLayers=None, contrastThreshold=None, edgeThreshold=None, sigma=None)参数的默认值分别为(nfeatures=0, nOctaveLayers=3, contrastThreshold=0.04, edgeThreshold=10, sigma=1.6)\n",
    "\t+ nfeatures表示最终保留的特征点数，由于未加参数时返回124，所以此值设置为超过124，就对特征提取结果没有影响了；小于124时，会有些角点未被检测到；\n",
    "\t+ contrastThreshold / nOctaveLayers = 2 * threshold / 255; 其中threshold为DoG空间里的像素阈值，由于在OpenCV里SIFT算子要求输入图像深度为CV_8U，因而threshold处于[0,255]闭区间内才有意义，因此 f = contrastThreshold / nOctaveLayers取值范围为[0,2]闭区间；按照默认值取值的话，f = 0.04 / 3 = 0.013.\n",
    "+ FAST角点检测的结果也比较精准，FAST的原理就是只着眼于局部，对于黑白棋盘这种灰度值全局分布和局部分布都比较统一的图像，即局部最优就是全局最优的情形，FAST的结果就不会受过于局部性而变差；另外FAST的运行速度最快。\n",
    "+ ORB角点检测的结果跟FAST一致，但因为多了几何不变性的处理，及特征描述子的计算过程，所以运行时间比FAST慢；但运行速度比SIFT快。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 对于自选风景图的运行结果\n",
    "+ Harris角点检测\n",
    "![Harris角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/scene_harris.png)\n",
    "\n",
    "+ SIFT角点检测\n",
    "![SIFT角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/scene_sift.png)\n",
    "\n",
    "+ FAST角点检测\n",
    "![FAST角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/scene_fast.png)\n",
    "\n",
    "+ ORB角点检测\n",
    "![ORB角点检测](https://gitee.com/coolhenry/cv-week-2/raw/master/out/blur/scene_orb.png)\n",
    "\n",
    "#### 对于自选风景图的结果分析\n",
    "+ Harris角点检测，检测的结果，取决于阈值的选取，容易出现响应过大的情况；\n",
    "+ SIFT角点检测，检测结果较好，使用默认阈值的适用性较好，而且具有尺度不变性；\n",
    "+ FAST角点检测，检测的结果，取决于图像本身灰度的分布情况，以及对特征点数量的设定；\n",
    "+ ORB角点检测，检测的结果较少。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
